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Abstract 

A C++ library, named ZKCM, has been developed for the purpose of multiprecision matrix calculations, which 
is based on the GNU MP and MPFR libraries. It is especially convenient for writing programs involving tensor- 
product operations, tracing-out operations, and singular-value decompositions. Its extension library, ZKCM_QC, 
for simulating quantum computing has been developed using the time-dependent matrix-product-state simulation 
method. This report gives a brief introduction to the libraries with sample programs. 
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1 Introduction 

Accuracy of simulation is often of serious concern when 
small differences in matrix elements result in physically 
important phenomena of one's interest. There are several 
programming libraries, e.g., Refs. [1, 2, 3, 4], useful for 
high-precision computing for this purpose. Among them, 
the library named ZKCM library [5], which I have been 
developing, is a C++ library for multiprecision complex- 
number matrix calculations. It provides several func- 
tionalities including singular value decomposition, ten- 
sor product calculation, and tracing-out operations. It is 
based on the GNU MP (GMP) [6] and MPFR [7] libraries, 
which are commonly included in recent distributions of 
UNIX-like systems. 

There is an extension library named ZKCM_QC. This li- 
brary is designed for simulating quantum computing [8, 9]. 
It uses a matrix product state (MPS) [10, 11] to represent 
a pure quantum state. The MPS method is recently one 
of the standard methods for simulation-physics software 
[12]. As for other methods effective for simulating quan- 
tum computing, see, e.g., Refs. [13, 14]. With ZKCM.QC, 
one may use quantum gates in U(2), U(4), and U(8) as ele- 
mentary gates. Indeed, in general, quantum gates in U(2) 
and U(4) are enough for universal quantum computing, 
but we regard quantum gates in U(8) also as elementary 
gates so as to reduce computational overheads. 

A simulation of quantum computing with MPS is known 
for its computational efficiency in case the Schmidt ranks 
are kept small during the simulation [11, 15]. Even for the 
case slightly large Schmidt ranks are involved, it is not as 
expensive as a simple simulation. This is obvious from the 
theory which is briefly explained in Sec. 3.1. 

This contribution is intended to provide a useful intro- 



duction for programming with the libraries. Section 2 de- 
scribes an example of simulating an NMR spectrum in a 
simple model using the ZKCM library. Section 3 shows 
an overview of the theory of the MPS method and an 
example of simulating a simple quantum circuit using the 
ZKCM_QC library. Effectiveness of the use of the libraries 
manifested by the examples are summarized in Sec. 4. 



2 ZKCM Library 

The ZKCM library is designed for general-purpose ma- 
trix calculations. This section concentrates on its main 
library. It consists of two major C++ classes: zkcm_class 
and zkcmjnatrix. The former class is a class of a complex 
number. Many operators like "+=" and functions like 
trigonometric functions are defined for the class. The lat- 
ter class is a class of a matrix. Standard operations and 
functions like matrix inversion are defined. In addition, 
the singular- value decomposition of a general matrix, the 
diagonalization of an Hermitian matrix, discrete Fourier 
transformation, etc., are defined for the class. A detailed 
document is placed in the "doc" directory of the pack- 
age of ZKCM. We will next look at a simple example to 
demonstrate the programming style using the library. 



2.1 Program example 

Here, a sample program "NMR_spectrum_simulation.cpp" 
found in the "samples" directory of the package of ZKCM 
is explained. This program generates a simulated FID 
spectrum of liquid-state NMR for the spin system consist- 
ing of a proton spin with precession frequency w\ = 400 
MHz (variable wl in the program) and a 13 C spin with 



precession frequency w 2 = 125 MHz (variable w2) at room 
temperature (300 K) (variable T). A J coupling constant 
J12 = 140 kHz (variable J12) is considered for the spins. 

The first line of the program is to include a header file 
of ZKCM: 

#include "zkcm.hpp" 

int main (int argc , char *argv[]) 

{ 

In the main function, the internal precision is set to 280 
bits for floating-point computation by 

zkcm_set_default_prec(280) ; 

In the subsequent lines, Pauli matrices /, X, Y, Z are gen- 
erated. For example, Y is generated as 

zkcm_matrix Y(2,2); 

Y. set (zkcm_class(0. 0,-1.0) ,0,1) ; 

Y.set(zkcm_class(0.0,1.0) ,1,0) ; 

Similarly, the Ygo pulse is generated as 

Yhpi . set (sqrt (zkcm_class (0 . 5) ) ,0 ,0) ; 
Yhpi . set (sqrt (zkcm_class (0.5)) ,0,1) ; 
Yhpi . set (-sqrt (zkcm_class (0 . 5) ) , 1 , 0) ; 
Yhpi . set (sqrt (zkcm_class (0 . 5) ) , 1 , 1) ; 

Other matrices are generated by similar lines. After this, 
values of constants and parameters are set. For example, 
the Boltzmann constant fee [J/K] is generated as 

zkcm_class kB("l .3806504e-23") ; 

(Several lines are omitted in this explanation.) The Hamil- 
tonian H in the type of zkcm_matrix is made as 

H = wl * tensorprod(Z , I) + w2 * tensorprod(I ,Z) 
+ J12 * tensorprod(Z.Z) ; 

This is used to generate a thermal state p: 

zkcm_matrix rho(4,4); 

rho = exp_H((-hplanck/kB/T) * H) ; 

rho /= trace (rho); 

Here, exp_H is a function to calculate the exponential of 
an Hermitian matrix and hplanck is the Planck constant 
(6.62606896 x 10~ 34 Js). The sampling time interval dt 
to record the value of < X > for the proton spin is set to 
0.145/iUi (any number sufficiently smaller than 1/2 might 
be fine instead of 0.145) by 

zkcm_class dt (zkcm_class ( "0 . 145")/wl) ; 

The number of data to record is then decided as 

int N = UNP2(1.0/dt/J12) ; 

Here, function UNP2 returns the integer upper nearest 
power of 2 for a given number. Now arrays to store data 
are prepared as row vectors. 

zkcm_matrix array (1, N) , array2(l,N); 

The following lines prepare the X, Y, and Fgo-pulse oper- 
ators acting only on the proton spin. 



zkcm_matrix XI (4, 4), Yl(4,4), Yhpil(4,4); 
XI = tensorprod(X, I) ; 
Yl = tensorprod(Y, I) ; 
Yhpil = tensorprod(Yhpi, I); 

To get an FID data, we firstly tilt the proton spin by the 
ideal pulse. 

rho = Yhpil * rho * adj oint (Yhpil) ; 

Now the data of time evolution of < X > of the proton 
spin under the Hamiltonian H is recorded for the time 
duration N x dt using 

array = rec_evol(rho, H, XI, dt , N) ; 

We now use a zero-padding for this "array" so as to en- 
hance the resolution. This will extend the array by N 
zeros. 

array2 = zero_padding (array, 2*N) ; 

To obtain a spectrum, the discrete Fourier transformation 
is applied. 

array2 = abs (DFT(array2) ) ; 

The "array2" is output to the file "example_zp.fid" as an 
FID data with df = 1/(2 x 2N x dt) as the frequency 
interval, in the Gnuplot style by 

GP_lD_print(array2, 1 . 0/dt/zkcm_class(2*N*2) , 
1, "example_zp.f id") ; 

At last, the function "main" ends with return 0;. The 
program is compiled and executed in a standard way. 1 

The result stored in "example_zp.fid" is visualized by 
Gnuplot as shown in Fig. 1. 
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Figure 1: Plot of the simulation data stored in "exam- 
ple_zp.fid" . See Software information for the environment 
where this simulation has been performed. 

It should be noted that we did not employ a high- 
temperature approximation. Under a high-temperature 
approximation, the first order deviation —f3H of 
exp(— f3H) (here, j3 = h/(k B T)) is considered as a devia- 
tion density matrix and calculations are performed using 



lr To make an executable file, the library flags typically "-lzkcm - 
1m -lmpfr -lgmp -lgmpxx" are probably required. As for ZKCM_QC, 
additionally "-lzkcm.qc" should be specified. 



the normalized deviation density matrix — if /const. This 
approximation is commonly used [17] but it cannot be 
used for simulations for low temperature. An advantage 
of using ZKCM for simulating NMR spectra is that the 
temperature can be chosen. This is possible because of 
high accuracy in computing the exponential of a Hamilto- 
nian. 

3 ZKCM_QC library 

The ZKCM_QC library is an extension of the ZKCM li- 
brary. It has several classes to handle tensor data useful for 
the time-dependent MPS simulation of a quantum circuit. 
Among the classes, the "mps" class and the "tensor2" class 
will be used by user-side programs. The former class con- 
ceals the complicated MPS simulation process and enables 
writing programs in a simple manner. The latter class is 
used to represent two-dimensional tensors which are often 
simply regarded as matrices. A quantum state during an 
MPS simulation is obtained as a (reduced) density matrix 
in the type of tensor2. For convenience, there is a function 
to convert a matrix in the type of tensor2 to the type of 
zkcm_matrix. 

More details of the classes are explained in the doc- 
ument placed at the "doc" directory of the ZKCM_QC 
package. 

We briefly overview the theory of the MPS simulation 
before introducing a program example since understand- 
ing the behavior of a library leads to a better programming 
in general. 

3.1 Brief overview of the theory of time- 
dependent MPS simulation 

Consider an n-qubit pure quantum state 
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with ^ i t i \ci a ...i n _ 1 2 = 1. If we keep this state as 
data as it is, updating the data for each time of uni- 
tary time evolution spends 0(2 2n ) floating-point opera- 
tions. To avoid such an exhaustive calculation, in the 
matrix-product-state method, the data is stored as a kind 
of compressed data. The state can be represented in the 
form 
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where we use tensors {Q S }™Z with parameters i s ,v s -i, v s 
(v-i and v n -i are excluded) and {V s }"Iq with parameter 
v s ; m s is a suitable number of values assigned to v s with 
which the state is represented precisely or well approxi- 
mated. This form is one of the forms of matrix product 



states (MPSs). The data are compressed to tensor ele- 
ments. We can see that neighboring tensors are correlated 
to each other; the data compression is owing to this struc- 
ture. 

Let us explain a little more details: Q s (i S7 v s -i,v s ) is 
a tensor with 2 x m s -i x m s elements; V s (v s ) is a tensor 
in which the Schmidt coefficients for the splitting between 
the sth site and the (s + l)th site (i.e., the positive square 
roots of non-zero eigenvalues of the reduced density oper- 
ator of qubits 0, . . . , s) are stored. This implies that, by 
using V s and eigenvectors |<&2"' s ) (|$* +1 "" n_1 )) of the re- 
duced density operator p°— s (p s +!— Q f qubits . . . , s 
(s + 1, . . . , n — 1), the state can also be written in the form 
of Schmidt decomposition 



s+l. 



(2) 



v s =0 



In an MPS simulation, very small coefficients and cor- 
responding eigenvectors are truncated out unlike a usual 
Schmidt decomposition. 

The advantage to use the MPS form is that we have 
only to handle a small number of tensors when we sim- 
ulate a time evolution under a single quantum gate. For 
example, when we apply a unitary operation £ U(4) act- 
ing on, say, qubits s and s + l, we have only to update the 
tensors Q s (i s ,v s -i,v s ), V s (v s ), and Q s+1 (i s+1 ,v s ,v s+ i). 
For the details of how tensors are updated, see Refs. 
[11, 16]. The simulation of a single quantum gate £ U(4) 
spends 0(m^ ax ) floating-point operations where m max is 
the largest value of m s among the sites s. (Usually, unitary 
operations £ U(2) and those £ U(4) are regarded as ele- 
mentary quantum gates.) A quantum circuit constructed 
by using at most g single-qubit and/or two-qubit quantum 
gates can be simulated within the cost of 0(gnmf nBX max ) 
floating-point operations, where n is the number of wires 
and TOmax,max is the largest value of m max over all time 
steps. 

The computational complexity may be slightly differ- 
ent for each software using MPS. In the ZKCM.QC li- 
brary, we have functions to apply quantum gates £ U(8) 
to three chosen qubits. Internally, three-qubit gates are 
handled as elementary gates. This makes the complex- 
ity a little larger. A simulation using the library spends 
0(gnm^ Illlx max ) floating-point operations, where g is the 
number of single-qubit, two-qubit, and/or three-qubit 
gates used for constructing a quantum circuit. 

The MPS simulation process, which is in fact often com- 
plicated, can be concealed by the use of ZKCM_QC. One 
may write a program for quantum circuit simulation in an 
intuitive manner. Here is a very simple example. 



3.2 Program example 

The following program is placed at the "samples" direc- 
tory of the ZKCM_QC package. It utilizes several matrices 
declared in the namespace "tensor2tools" (see the docu- 
ment for details on this namespace). 

#include "zkcm_qc .hpp" 



int main (int argc, char *argv[]) 
{ 

//Use the 256-bit float for internal computation. 
zkcm_set_def ault_prec(256) ; 

//Num. of digits for each output is set to 8. 
zkcm_set_output_nd(8) ; 

//First, we make an MPS representing I 000> . 
mps M(3) ; 

std::cout « "The inital state is " « std::endl; 
//Print the reduced density operator of the block 
//of qubits from to 2, namely, 0,1,2, using the 
//binary number representation for basis vectors. 
tensor2tools: : showb(M.RD0_block(0, 2) ) ; 

std::cout « "Now we apply H to the 0th qubit . " 

« std: :endl; 
M.applylI(tensor2tools: :Hadamard, 0) ; 



ZKCM.QC version 0.0.1 on the Fedora 15 64-bit oper- 
ating system with GMP version 4.3.2 and MPFR version 
3.0.0 have been used for this report. 
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